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Abstract 

Several models of Quantum Gravity predict Lorentz Symmetry breaking at en- 
ergy scales approaching the Planck scale (~ 10 19 GeV). With present photon 
data from the observations of distant astrophysical sources, it is possible to 
constrain the Lorentz Symmetry breaking linear term in the standard photon 
dispersion relations. Gamma-ray Bursts (GRB) and flaring Active Galactic 
Nuclei (AGN) are complementary to each other for this purpose, since they 
are observed at different distances in different energy ranges and with differ- 
ent levels of variability. Following a previous publication of the High Energy 
Stereoscopic System (H.E.S.S.) collaboration a more sensitive event-by-event 
method consisting of a likelihood fit is applied to PKS 2155-304 flare data of 
MJD 53944 (July 28, 2006) as used in the previous publication. The previous 
limit on the linear term is improved by a factor of ^3 up to Mq G > 2.1 x 10 18 
GeV and is currently the best result obtained with blazars. The sensitivity to 
the quadratic term is lower and provides a limit of Mq G > 6.4 x 10 10 GeV, which 
is the best value obtained so far with an AGN and similar to the best limits 
obtained with GRB. 

Keywords: active galaxies, PKS 2155-304, H.E.S.S., quantum gravity, 
Lorentz invariance breaking 



1. Introduction 

1.1. Lorentz Invariance and Quantum Gravity 

A development of new space-borne and ground-based experiments in the 
last decade, covering a large domain of gamma-ray astronomy, opened a new 
window for fundamental physics tests. In particular, the confrontation of the 
results from astrophysical sources with challenging predictions of various models 
in the frame of Quantum Gravity theories became of general interest [2J . The 
theory of Quantum Gravity (see e.g. Q), which is still incomplete, would give a 
unified picture based on both Quantum Mechanics and General Relativity, thus 
leading to a common description of the four fundamental forces. 

Quantum Gravity effects in the framework of String Theory where the 
gravitation is considered as a gauge interaction, result from a graviton-like ex- 
change in a classical space-time. In most String Theory models involving large 
extra-dimensions, these effects would take place at the Planck scale, thus not 
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leading to "spontaneous" Lorentz Symmetry breaking, as may happen in mod- 
els with "foamy" structure @, 0] of quantum space-time. In this second class 
of models, photons propagate in a vacuum which may exhibit a non-negligible 
refractive index due to its foamy structure on a characteristic scale approaching 
the Planck length or equivalently the Planck energy Ep = 1.22 x 10 19 GeV. 
This implies a group velocity of light changing as a function of energy of the 
photons, in analogy to the dispersion effects in any theoretical description of 
plasma media. 

On the other hand, in models based on General Relativity with Loop Quan- 
tum Gravity [H 0] which postulate discrete (cellular) space-time in the Planck- 
ian regime, the fluctuations would introduce perturbations to the propagation 
of photons. The light wave going through the discrete space matrix would feel 
an induced perturbation which increases with decreasing wavelength or equiva- 
lently with increasing energy of photons. In consequence, photons with different 
wavelengths propagate with different velocity. 

As a result, one may expect a spontaneous violation of Lorentz Symmetry 
at high energies to be the generic signature of the Quantum Gravity. 

In four dimensions, the Quantum Gravity scale is presumed to be close to the 
Planck mass and the standard photon dispersion relations up-to second order 
corrections in energy can be written as: 

cV = E 2 (1 ± i(E/M) ± ((E/M) 2 ± ...) , (1) 

where £ and C are positive parameters. 

This implies an energy-dependent speed-of-light with 

« = c(l-f(£7/M)) (2) 

considering only the linear term in the expansion. Alternatively assuming the 
linear correction equal to zero, the quadratic term leads to: 

v = c(l-((E/M) 2 ). (3) 

High-energy photons could propagate either slower (the sub-luminal case) 
or faster (the super-luminal case) than low-energy photons. In this study, a 
causality conserving scheme was adopted following to the theoretical arguments 
of the non-critical string models, developed in Q. The sign minus in equations 
[5] and [3] is required to avoid Cherenkov radiation in vacuum [9( . As results of 
this analysis and as in case of most of other quoted results in Table [TJ only 
limits on Quantum Gravity scale for the sub-luminal photons will be given. 

As suggested by Amelino-Camelia [5|, the tiny effects due to Quantum Grav- 
ity can add up to measurable time delays for photons from cosmological sources. 
The energy dispersion is best observed in sources that show fast flux variability, 
that are at cosmological distances and are observed over a wide energy range. 
This is the case of Gamma Ray Bursts (GRB) and Very High Energy (VHE) 
flares of Active Galactic Nuclei (AGN). Both types of sources are the preferred 
targets of these time-of-fiight studies, which provide the least model-dependent 
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tests of the Lorentz Symmetry. The case of pulsed emission by Pulsars has also 
been considered, and provided valuable results flo| ] . 

It should be underlined that studying both GRB and AGN is of fundamental 
interest. GRB are observed with satellites at very large distances (up to z ~ 8) 
but with very limited statistics above a few tens of GeV. On the contrary, 
AGN flares can be detected by ground based instruments with large statistics 
up to a few tens of TeV. Due to the absorption of the high energy photons by 
Extragalactic Background Light (EBL), TeV observations are limited to sources 
with low redshifts. Thus, GRB and AGN are complementary since they allow 
to test different energy and redshift ranges. 

It should be noted that time-of-flight measurements are subject to a bias 
related to a potential dispersion introduced from the intrinsic source effects, 
which could cancel-out or enhance the dispersion due to modifications to the 
speed of light. 

When considering sources at cosmological distances, the analysis of time lags 
as a function of redshift requires a correction due to the expansion of the universe 



which depends on the cosmological model. Following the analysis of the 



BATSE data and recently of the HETE-2 and Swift GRB data [jl O, Q [ljj 
and within a framework of the Standard Cosmological Model [16| with a flat 
expanding universe and a cosmological constant, the difference in the arrival 
time of two photons emitted at the same time is given by 

At £ f z (l + z')dz' 

I (4) 



(5) 



AE E P H o y o ^ m (l + 2 ') 3 + A 
for a linear effect, and by 

At 3C [ z {l + z'fdz 1 



AE* 2E2H y ^(I + z'^ + Oa 



for a quadratic effect. In the present study the cosmological parameters were 
set to Q m = 0.3, ft A = 0.7 and H = 2.3 x KT 18 s -1 . AE and AE 2 represent 
linear and quadratic energy ranges respectively. 

1.2. Present results 

The most important results obtained so far with transient astrophysical 
sources are given in Table Q1 The most robust limits to date have been ob- 



tained by Ellis et al. 13|, [l4| and Bolmont et al. 15j, with the use of several 
GRB at different redshifts, resulting in limits of the order of Mq G > 10 16 GeV. 
In case of a significant detection, using several sources at different redshifts al- 
lows in principle to take into account the possible time-lags originating from 
source effects. 



Fermi is now in operation and has provided the best limits so far 17J, LL8J . 
However, for the moment the results are obtained with no study of the variation 
of the time-lag with redshift and with limited statistics in the GeV range. The 
recent studies with flares of Mrk 501 19|, 2(| and the former result of H.E.S.S. 
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Table 1: A selection of limits obtained with various instruments and methods for both GRB and AGN 



Source(s) 



Experiment 



Method 



Results (95% CL limits) 
> 1.8xl0 17 GeV 



Reference Note 



Mqg 

Kg 

Kg 

Kg 

Kg 

Kg 

Kg 

Kg 
M] 



tell 


a b 


[17] 




[18] 


c 
c 


[121 


a 


[15] 


d 


[22] 


e 


[13, M 


f g 

1 


m 


a h 
j 


[19] 




[20] 




[1] 





GRB 021206 
GRB 080916C 



RHESSI 

Fermi GBM + LAT 



GRB 090510 Fermi GBM + LAT 



9 GRBs 

15 GRBs 
17 GRBs 
35 GRBs 

Mrk 421 
Mrk 501 



BATSE + OSSE 

HETE-2 
INTEGRAL 
BATSE + HETE-2 
+ Swift 
Whipple 
MAGIC 



PKS 2155-304 H.E.S.S. 



Fit + mean arrival time in a spike 
associating a 13 GeV photon with the 
trigger time 

associating a 31 GeV photon with the 
start of any observed emission, DisCan 
wavelets 

wavelets 

likelihood 

wavelets 

likelihood 
ECF 

likelihood 

MCCF 

wavelets 
likelihood 



M: 



QG 
QG 



> 1.3xl0 18 GeV 

> 0.8xl0 10 GeV 

> 1.5xl0 19 GeV 

> 3.0xl0 10 GeV 

> 0.7xl0 16 GeV 

> 2.9xl0 6 GeV 

> 0.4xl0 16 GeV 

> 3.2xl0 n GeV 

> 1.4xl0 16 GeV 



Kg 
M ; qg 

Kg 

Kg 
Kg 
Kg 
Kg 
Kg 
Kg 
Kg 



> 0.4x 

> 0.2x 

> 0.3x 

> 0.3x 

> 5.7x 

> 7.2x 

> 1.4x 

> 5.2x 

> 2.1x 

> 6.4x 



10 17 GeV 

10 18 GeV 
10 11 GeV 
10 18 GeV 
10 10 GeV 
10 17 GeV 
10 9 GeV 

10 17 GeV 

10 18 GeV 
10 10 GeV 



"Limit obtained not taking into account the factor (1 + z) in the integral of equations 0] and [5] b The pseudo-redshift estimator 
32j was used. This estimator can be wrong by a factor of 2. c Only the most conservative limits are given here. The best 
limits given in 18 1 which can still be considered as secure are M l QG > 10 E P and Mq G > 8.8 xlO 10 GeV. d Photon tagged data 
was used. e The pseudo-redshift estimator [32| was used for 6 GRB out of 11. •'For HETE-2, fixed energy bands were used. 
9 The limits of Ellis et al. [l3| were updated in [l4| taking into account the factor (1 + z) in the integral of equation |U Only 
the limit obtained for a linear correction is given. h A likelihood procedure was used, but not on an event- by- event basis. l This 
work. 



with PKS 2155-304 [l| provide results setting a lower limit on the linear Quan- 
tum Gravity scale at a few 10 17 GeV for each source. 

Table [T] also shows the different analysis techniques used to measure the 
energy-depend time-lags. The analysis methods can be divided in two categories, 
analyses applied on binned data on the one hand and unbinned analyses on the 
other hand. 

Binned analyses are applied on binned light curves and consist of finding 
the position of local extrema in two different energy bands and deducing the 
value of the lag. A simple fit [2l[ or a more sophisticated technique such as 
Continuous Wavelet Transform (CWT) can be used to localize extrema. CWT 
was used for GRB [13, El 111 111 and PKS 2155-304 Q. The Modified Cross 
Correlation Function (MCCF) was also used for PKS 2155-304 to determine 
directly the lag from two light curves in two different energy bands [if. 

Unbinned analyses have also been used for both GRB and AGN. Likelihood 
fits require a parameterization of the light curves and this parameterization can 
be taken from a binned 20] or unbinned fit 22] to the light curve. Another 
method, the Energy Cost Function (ECF), was recently used by Albert et al. 
(l9j |. In this method one introduces a time lag r for each photon and adds up 
their energies if they fall into a given time interval. The maximum of ECF(r) 
gives the value of the measured lag. The recently developed method called 
DisCan [23| was used for the analysis of Fermi data [l8j] . As discussed in (23[ , 
this method applies successfully to the analysis of very sharp peaks as can be 
found in GRB light curves. 

It is important to point out here that the use of at least two methods is 
recommended since different methods can probe different aspects of the light 
curve. The present analysis completes previously published studies [H where 
two different methods were employed. 

As no significant time-lag was detected so far, careful error calibration studies 
using Monte Carlo simulations are mandatory for the extraction of the limits, 
as will be further discussed in the following. 

1.3. Overview of the paper 

In this paper, a new analysis of PKS 2155-304 data is presented, using a 
likelihood fit to determine both linear and quadratic corrections to the dispersion 
relations. In section^ the method is described. The data in use as well as the 
selections applied are detailed in section [3j Then, in section [4] the procedure 
used to calibrate the method and to evaluate systematics is presented. The 
results are given in section [5] and finally discussed in section [6] 

2. Method description 

To study the correlations between the arrival time and the energy of the 
photons, a method relying on a probability density function (p.d.f.) on an 
event-by-event basis is used. Following Martinez and Errando [20(, the p.d.f. is 
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defined by: 



P(t,E)=N A{E s )A{Es)G(E-Es,*{Es))F s {t-TiEs)dEs, (6) 
Jo 

where A(Es) is the emission spectrum, G(E — E$, a(Es)) is the smearing func- 
tion in energy (a usual point spread function for an energy resolution of 10%), 
A(Es) is the acceptance of the H.E.S.S. array, Fs is a parameterization of the 
emission time distribution (or light curve) and N is a normalization factor, the 
time shift was re-parameterized as At = —t\E and At = —r q E 2 for linear and 
quadratic effects respectively, where 77 and r q are expressed in s TeV -1 and 
s TeV -2 . The value of N was computed in view of unbiased T-parameter es- 
timation. It guarantees the correct normalization of the p.d.f. and takes into 
account its model dependence. Here, integration over measured domain of en- 
ergy and time has been performed. 

The function Fs is a function of time at the source tg = t — tEs- The 
probability P(t, E) depends only on the relative delay of photons for different 
energies. As the light curve at the emission time is not known, the best estimate 
of the time structure at the source can be taken from the time distribution of 
the lowest energy photons. In the propagation models, the change of the QG 
scale induces basically a sliding effect in time of the light curve. This point has 
been checked to be valid by the Monte Carlo simulations described below. 

An important assumption of this study should be underlined: the intrinsic 
light curves have been considered to have the same shape in different energy 
bands. As no lag is found, this can be partially justified by comparing the light 
curve at low energy (Fig. ^ and in the full energy range of the instrument 



where the statistics is five times higher [24[ . In both cases the light curve shows 



a five peak structure with maxima located roughly at the same times. 

The likelihood function is given by the product of probability density func- 
tions over all photons: 

L = l[P t (t,E). (7) 

i 

The maximum of L provides the time-lag parameter 77 and r q for linear and 
quadratic models respectively. 

The likelihood method is in principle very sensitive and allows a lag to be 
measured even if the number of selected photons is limited, which is not the 
case for wavelet or MCCF analyses. As a consequence of this high statistical 
sensitivity, it would be possible to probe a light curve locally, introducing a 
selection in time. The main aspect of the method to be taken into account is a 
need for a very precise parameterization of the photon distributions. 



3. Data sample used in the analysis 

PKS 2155-304 is a high-peaked BL Lac object located at redshift z = 0.116. 
On July 28, 2006 (MJD 53944), an extreme flare of this source was observed by 
H.E.S.S. [24|. The full data sample was obtained using a combination of two 
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Tabic 2: Selections applied to the data 

Selection Number of selected events 



Total sample 




8148 (100 %) 


(1) = Time in 


0-4000 s 


7713 (95 %) 


(2) = (1) and 


9 2 < 0.005 deg. 2 


5974 (73 %) 


(2) and E m in 


0.25-4.0 TeV 


3526 (43 %) 


(2) and E m in 


0.25-0.28 TeV 


561 (7 %) 



analysis methods (25[: the standard Hillas analysis [26| and the Model analysis 
27| • This technique greatly improves the signal to background ratio and results 
in 8148 on-source events recorded in ~85 minutes (three data-taking runs) with 
an energy threshold of ~120 GeV. The signal-to-background ratio in the flare 
data exceeds 300. The effect of the background contribution will be investigated 
in the next section. The average zenith angle of selected events is Z ~ 11°. 

Table [2] summarizes the different cuts applied to the data. In this study, 
only the first 4000 s of data (7713 events) are considered, where the flux and 
variability are the highest. To further improve the signal-to-background ratio, a 
cut 9 2 < 0.005 deg. 2 is applied. This cut reduces the number of photons to 5974 
events but allows a very good fit of the light curve and of the energy spectrum. 
A total of 3526 events pass the cut on 9 2 in the first 4000 s of data in the range 
0.25-4.0 TeV. All events in the data sample fulfilling the selection were used in 
the analysis. 

As shown in Figure[TJ the reconstructed spectrum is compatible with a power 
law of index T = 3.46 in the range of 0.25-4.0 TeV (with x 2 / dof = 16.7/23). 
This result is compatible with the high energy part of the broken power law 
spectrum obtained for the same flare by Aharonian et al. [24| with a different 
event selection and energy reconstruction method, and also with the spectrum 
obtained for the whole 2005-2007 period [28j]. The dependence of the spectral 
index as a function of time has been studied extensively in [24| and no significant 
variation was found during the flare. 

Figure [2] shows the light curve in the range 0.25-0.28 TeV with a binning 
of 61 s. The shape of this light curve will be used later as a template for the 
time-lag determination (31]). To estimate the main parameters of the light curve 
(widths and asymmetry of the pulses, distance between consecutive spikes), a 
parameterization was performed with a sum of asymmetric Gaussian functions. 
For this, the number of spikes was first determined using a peak-finding proce- 
dure based on a Markov chain algorithm (29| . Five pulses were found and their 
positions were used as initial parameters for the light curve fit. 

The fit function F is a sum of asymmetric Gaussian curves, defined as follows: 

71 

F n {t)=^f i (t,A i ,i H ,a i ,p i ), (8) 
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Figure 1: Spectrum of PKS 2155-304 flare of MJD 53944, taking into account the cut on the 
off-axis angle mentioned in the text. Points are fitted with a power law E~ T . 
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Figure 2: Light curve of PKS 2155-304 flare of MJD 53944 in the range 0.25-0.28 TeV, with 
a bin width of 61 s and taking into account the cut on 9 2 described in the text. Black points 
show the positions of the extrema as determined by the peak-finding procedure of Morhac et 
al. |29l . The light curve is fitted with a function F$(t) + B where B is a constant and F§ is 
defined by equation [8] The zero of the time axis corresponds to MJD 53944.024. 
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Table 3: Parameters of the function i*g(t) + B fitted on PKS 2155-304 light curve in the range 
0.25-0.28 TeV 



Parameter 


Value 


Error 


Ai 


27.1 


1.4 


Mi 


430.4 


35.8 


ai 


399.9 


25.4 


Pi 


352.9 


75.7 


A 2 


26.0 


0.8 




1456.1 


16.0 




342.9 


92.8 




193.3 


44.9 




zo.v 


U. 1 


M3 


2113.1 


19.2 




276.1 


46.2 


& 


296.8 


29.3 


A 4 


20.0 


0.7 


fi 4 


2639.7 


77.0 


04 


130.3 


34.3 




899.8 


23.7 


M 


9.4 


1.3 




3302.3 


19.9 


«5 


87.1 


36.5 


A> 


999.9 


80.4 


B 


-14.1 


0.8 



where n is the number of spikes, and where 

! («-«*)" 
y = A e 2^ , if t < /it 
_ (t -,) 2 ( 9 ) 
y = A e 2 f' 2 , if t > fi, 

In this definition, A and /1 are the normalization factor and the position of 
the pulse and a and ft are the left-hand and right-hand widths of the pulse 
respectively. The asymmetry of the curve can be quantified by the parameter 
v — ft I a. When v = 1, / becomes a Gaussian function with a width a — ft = oq. 

The parameterization was carried out with a function F, 5 plus a constant term 
B. Values of the parameters are given in Table [3j leading to x 2 /dof = 48.8/38. 
This fit gives the following information about the light curve: 

• width of the pulses in the range 250-600 s, 

• time difference between consecutive spikes in the range 500-1000 s, and 

• ratio v in the range 0.6-7. 

These intervals were investigated for the calibration procedure described in 
the next section. 
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4. Calibration of the method 

A calibration procedure was developed using a toy Monte-Carlo program in 
order to reflect as closely as possible the real data analysis and to provide the 
best strategy to minimize the systematic errors on the estimated time lag. 

As discussed in the previous section, the light curve of PKS 2155-304 is 
well fitted by a superposition of five asymmetric Gaussian spikes. The energy 
spectrum follows a power law distribution. These parameters have been used 
as initial conditions for the simulation. The zenith angle was considered to be 
constant at Z = 10°. 

A simulation procedure comprises the following steps: 

• the template F$ is chosen to be either a Gaussian function or derived from 
the measured light curve in the range 0.25-0.28 TeV; 

• photons are generated from a parameterizations F$ and A in the range 
0.25-4 TeV; 

• for each photon, a lag At is introduced depending on the energy; 

• a weight is associated to each photon taking into account the acceptance 
of the detector during data-taking in 2006; 

• the energy is smeared taking into account the H.E.S.S. energy resolution 
(A£/£~10%); 

• when a simple Gaussian light curve is generated, histograms obtained for 
the light curve and the spectrum are fitted. As a result, the functions 
used for the fits reproduce well those used for the generation. This point 
was checked to estimate the effects of the statistical smearing, which were 
found to be negligible; 

• finally, the likelihood is computed and the minimum of — 2Aln(L) is 
searched for, first using a Gaussian template parameterization for the 
detailed studies or the PKS 2155-304 template. As shown later, the dis- 
tribution of the minima has always been found to be compatible with a 
Gaussian distribution (with a mean f and a width er T ). 

For each simulation run, these steps are repeated 500 times with ~2000 
photons after selections in each realization. The injected lag (r ) range was 
chosen between —100 and 100 sTeV -1 (respectively sTeV -2 ) with steps of 
10 sTeV" 1 (sTeV- 2 ). 

The results of the simulations can be summarized studying the variation of 
the reconstructed lag r as a function of the injected lag. Figure [3] (left) shows 
such a plot for simulation settings which will be described later (Gaussian spike). 

The shaded area corresponds to a 68% CL range of the reconstructed lag. 
A fit within ±ov yields a slope parameter value close to unity. 
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Figure 3: Left: reconstructed lag as a function of the injected lag for a Gaussian light curve 
in the modelling mode (68% CL range). The plot in the bottom right corner shows the 
distribution of the minima for the 500 realizations with a generated lag of sTeV - . The 
dashed line shows the linear response function. Right: parameter a T as a function of the width 
of the Gaussian pulse. The axis at the top gives the full width at half maximum (FWHM) of 
the pulse. 

A small systematic shift of the slope is observed; however, the bias on the 
reconstructed time lag is much smaller than the statistical uncertainty a r for 
time lags between —50 sTeV -1 and 50 sTeV -1 . 

In the following sub-sections, different features of the error calibration proce- 
dure are detailed. First, the influence of the pulse shape (width and asymmetry) 
and the superposition of pulses is investigated. Then, the dependence on the 
spectral index is evaluated. In the end, the effect of background events is stud- 
ied. The simulation of light curves with five pulses provide the value of the 
calibrated statistical error which is used in £[5]to compute the limits on Eqg- 

4--1- Pulse shape and superposition of pulses 

Figure [3] (left) shows the reconstructed lag as a function of the injected 
lag for a Gaussian pulse of a mean width ctq of 300 s for a linear model. The 
calibration slope is 0.97±0.02, showing the likelihood fit provides efficient means 
for the reconstruction of a lag. No deviation from linearity is observed in this 
case. The same figure shows an example of the distribution of the minima 
of — 2Aln(L) fitted with a Gaussian curve, for an injected lag of sTeV -1 . 
This distribution is fitted by a Gaussian curve with f — 0.6 ± 0.2 sTeV -1 and 
a T = 5.1 ±0.2 sTeV" 1 . 

For a given value of the pulse width ctq > the spread a T is very stable. On the 
other hand, as shown in Figure [3] (right), a T increases rapidly with the width of 
the pulse. For a single pulse width of 600 s, a T rises to ~ 13 sTeV -1 . 

When assuming asymmetry in the pulse shape, the calibration slope remains 
very stable and close to unity. The maximal asymmetry factor observed for 
PKS 2155-304 is v = 7, leading to a T « 10 sTeV^ 1 for cr of 300 s. 
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Figure 4: Reconstructed lag (68% CL range) as a function of the injected lag for a light curve 
+ B obtained from real data and for a linear (left) and a quadratic (right) effect. The 
dashed line shows the linear response function. 

The superposition of spikes requires detailed studies concerning the evolu- 
tion of the calibration slope and the error in case of complex light curves such 
as the one measured with PKS 2155-304 flare. In order to check the effect of 
spike superposition, light curves were generated from the sum of several asym- 
metric Gaussian curves. Various configurations were tested and the resulting 
calibration curves were compared to those obtained when injecting the real 
PKS 2155-304 light curve. 

The plots of Figure U] were obtained with light curves generated from real 
data in the range 0.25-0.28 TeV for the linear and quadratic modelling. These 
simulations were used to investigate the effect of superimposing two to five pulses 
with realistic parameters leading to similar calibrated error values. 

In the quadratic case, the slope is somewhat smaller than 1. The spread a T 
is smaller than in the linear case close to To = but increases at larger values 
of |to[ where the calibration curve becomes non-linear. 

In the linear case, a systematic shift of the calibration values of the order of 
l/2<7 r is observed. This effect will be included in the overall calculation of the 
systematic errors (see Table HJ). 

The error values obtained from Figure 2] were used in the following to derive 
the limits for linear and quadratic corrections. The error <r T was found to depend 
mainly on individual pulse widths rather than on the addition of the widths of 
the five pulses, leading to 10.9 s TeV -1 for the linear and 6.3 s TeV -2 for the 
quadratic case respectively. Anticipating the results from data, these values 
were taken in the vicinity of the injected time-lag tq = 0. 

In addition, other functional forms were used when fitting template light 
curve: Norris parameterization as also used in [2~ij | and Landau function. No 
detectable effect was found when deriving the results. Moreover, the effect of 
the different functional forms was evaluated with MC simulations. 
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Figure 5: Parameter cr T as a function of the spectral index. 
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4-2. Dependence on spectral index 

Varying the spectral index of the generated light curves from 1 to 4.5 re- 
sulted in no significant change in the calibration slope (below 5%) for a single 
Gaussian pulse. The systematic error due to the uncertainty on the spectral 
index determination can then be considered as negligible. In contrast, the er- 
ror value is a strong function of the spectral index as also shown in Figure [SJ 
A harder spectrum gives a higher precision, due to the fact that the photons 
are spread more uniformly between high and low energies. This enhances the 
potential of studies with closer AGN such as Mrk 421 and Mrk 501. 

4-3. Effect of background contribution 

Even though the signal to background ratio is very high for the studied 
data sample, it is possible that a small contribution of mis-identified photons 
introduces a small bias in the measurement of the time-lag. To study the impact 
of this effect, light curves were simulated with a fraction of particles with no 
energy-dependant time-lag as is expected for charged cosmic-rays. A realistic 
power law spectrum with an index of 2.7 was used to generate these particles. 
When 1% of the particles are considered as protons, no significant change in 
either slope or parameter a T is observed. 

4-4- Summary of calibration studies 

In this section, the most important features of the calibrations with the 
likelihood method are summarized. 

As the emission light curve at the source is not known, the function Fs was 
derived from real data in the range 0.25-0.28 TeV. In this case, it is shown that 
the likelihood fit can reconstruct the injected lag extremely well. In the future, 
if a satisfying emission model is proposed, the likelihood approach can be very 
useful to test propagation and emission models as well. 

The error on the measured lag (a T ) depends strongly on the width of indi- 
vidual pulses. This result is very important and allows competitive results to be 
obtained with the PKS 2155-304 light curve. On the other hand, the asymmetry 
of the pulses and the superposition of several pulses has been shown to have less 
important effect on the results. The parameter a T also depends on the spectral 
index, which is another source of uncertainty. However, for PKS 2155-304 flare 
in 2006, the spectral index can be determined with a very good precision. 

The systematic uncertainties in the method employed were quantified by 
varying selections and cuts and treating the light curve shape parameterization 
in different ways. The various contributions are detailed in Table 01 The main 
uncertainty was found to be related to event selections and to the light curve 
parameterizations . 

In order to get a precise estimation of the systematic uncertainty related to 
the Fs parameterization, the covariance matrix of the fit presented in Figure [2] 
was studied. The values of the fit parameters were varied following to the 
gaussian distributions with parameters derived from the covariance matrix after 
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Table 4: Systematic uncertainties of the event- by- event likelihood 



Estimated error Change in estimated r\ (s TcV l ) Change in estimated r q (s TeV z ) 



Selection cuts 




< 5 


< 5 


Background contribution 


1% 


< 1 


< 1 


Acceptance factors 


2% 


< 1 


< 1 


Energy resolution 


1% 


< 1 


< 1 


Energy calibration 


10% 


< 2 


< 2 


Spectral index 


1% 


< 1 


< 1 


Calibration systematics (constant, shift) 


10% 


< 5 


< 1 


F s {t) parameterization 




w7 


«3 


Total 




< 10.3 


< 6.6 



Number of passing events = 2462 

— Minimum at x, n = -5.5 + 7.1 s/TeV 

< 




20 40 60 

x, (s TeV 1 ) 



Number of passing events = 2462 
Minimum at x nn = 1.7 + 3.5 s TeV" 2 




Us TeV 2 ) 



Figure 6: — 2Aln(L) as a function of r when the measured light curve is fitted in 0.25-0.28 
TeV and the likelihood is computed in 0.3—4.0 TeV for a linear (left) and quadratic (right) 
correction. Points are fitted with a third-degree polynomial with a minimum at t;o = — 5.5±7.f 
sTeV -1 and r q o = f.7 ± 3.5 sTeV - 2 . The errors on these values are obtained requesting 
-2Aln(L) = f. 



its diagonalization. The values in Tabled were obtained with this procedure for 
the linear and the quadratic case separately. 

The overall systematic error is estimated to be < 10.3 s TeV -1 and < 6.6 s TeV 
for the linear and quadratic effects respectively. These values will be combined 
with the statistical calibrated error given in 34. II when determining the Quantum 
Gravity scales. 

5. Results with PKS 2155-304 flare of MJD 53944 

The final result was obtained with a parameterization of the light curve 
in the range 0.25-0.28 TeV (El = 0.26 TeV, see light curve in Figure [2]) and 
a likelihood computed in 0.3-4.0 TeV (El = 0.49 TeV). The choice of these 
ranges was dictated by the balance found between lever-arm in energy and the 
statistics used in the likelihood fit. These choices were studied and optimized 
with the simulations described above. 

Figure E] shows the — 2Aln(L) curve as a function of the injected r for the 
linear and quadratic cases, leading to the value for the minimum of 

r u; = -5.5 ± 10.9 (stat) ± 10.3 (sys) sTeV" 1 

for the linear and 

r' 0q = 1.7 ± 6.3 (stat) ± 6.6 (sys) s TeV" 2 
for the quadratic effect. 
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These results imply that there is no significant time-lag measured. Conse- 
quently, the 95% CL limits on the Quantum Gravity energy scales are derived: 

Mq G > 2.1 x 10 18 GeV (£ < 5.7) 

and 

Mq G > 6.4 x 10 10 GeV (C < 3.6 x 10 16 ). 
The more conservative 99% CL limits are 

Kg 99% > 1-3 x 10 18 GeV (£ < 9.2) 

and 

M QG 99% > 53 x 1Ql ° GeV (C < 5 - 2 x lo16 )- 
These results were found to be very stable when considering different energy 
domains and time selections. The stability studies of the estimated systematic 
effects allows the quoted values to be considered as conservative. 

The limit on the linear term, while lower by a factor of ten in comparison 
with the one obtained by Fermi with GRB 090510 [3, shows that with AGN 
too, the conclusion about the models predicting a linear effect with energy 0] 
tends to be confirmed. As far as the limit on the quadratic term is concerned, 
it is the best one obtained so far with AGN and GRB, even if it still remains 
far from the Planck scale value. 

As stated in the introductory part, the presented paper focused on models 
proposed in [liL [l3| where predictions of linear Lorentz Invariance Violation 
effects were expected and the causality condition was conserved. Moreover, it 
may be noticed that because of negligible negative values found both for the 
linear and quadratic terms, the super-luminal case limits will differ only slightly 
from those provided for the sub-luminal case. 

6. Summary and prospects 

The limits obtained in this study are the most constraining ever obtained 
with an AGN. They are almost ten times higher than those obtained by Martinez 
and Errando [2(| with the same method for the Mrk 501 flare recorded by 
MAGIC, and consistent with their rough estimate for the PKS 2155-304 flare 
observed by H.E.S.S.. They are also a factor of ^3 higher than the previous 
H.E.S.S. result [l|. The increase in sensitivity as compared to the MAGIC result 
is mainly due to excellent parameters of the data taken on July 28, 2006: low 
zenith angle (~ 10°) which leads to a low energy threshold and high statistics 
with high variability. On the other hand, the steepness of the PKS 2155-304 
energy spectrum was one of the penalizing factors in this study. 

Given the complexity of the PKS 2155-304 light curve, a detailed calibra- 
tion study of the method was necessary. For this purpose, a toy Monte-Carlo 
program was developed. The simulations allowed the influence of different pa- 
rameters on the likelihood performance to be investigated, such as the shape of 
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the light curve or spike superposition. In particular, a very interesting conclu- 
sion for a light curve with multiple spikes is that the error on the lag depends 
mainly on the width of the individual pulses. This feature of the likelihood fit 
allows a smaller error on the lag measurement to be obtained, yielding higher 
limits on Mqq- 

It should still be noticed that the superposition of different spikes leads 
to a decrease in sensitivity due to a limited number of "useful events" which 
contribute to the likelihood fit. Despite this loss in statistical power, the method 
allows robust and precise results on the fitted time-lag to be obtained. 

Another important conclusion about the method is the possibility to obtain 
robust results even with low photon statistics. The use of different time selec- 
tions showed that a large decrease in the number of selected events entering the 
fit had only a small influence on the error value of the time-lag measurement. 

The main difficulty of time of flight studies with photons from AGN and 
GRB is that the origin of the measured lag is unknown. If the measured time- 
lag is not related to Quantum Gravity induced propagation effect, it can be due 
to an emission effect at the source. The measured lag can also be a superposition 
of source and propagation effects. As no lag was found in the studied sample, 
either the source effects compensate exactly the propagation effects ("conspir- 
acy scenario"), or there are no source or propagation induced effects detected 
within present experimental sensitivity. In this study the second hypothesis was 
retained, which leads to a conclusion that any measured lag would be entirely 
due to Lorentz Invariance Violation effects. 

In future, modelling of the emission light curve could bring new insight in the 
propagation and emission time studies. Population studies of AGN examining 
the effect of the distance on the lag measurement could allow this assumption 
to be avoided but would require that intrinsic effects are similar for different 
AGN. Better results on propagation effects will then come in parallel with a 
better understanding of the emission models. 

More and more results on Lorentz symmetry breaking are published which 
give limits for the linear correction that are very close to the Planck energy 
scale. In particular, the latest results from Fermi observations of distant GRB 
provide exclusion limits for a large class of linear models above the Planck energy 
[17l . [l8j . However, as pointed out by Ellis et al. [3(J, 'extraordinary claims 
require extraordinary evidence', so more studies will be needed in the future to 
give more robust results and to be able to definitively reject or validate proposed 
models. 

The time-lag measurements with photons are far from being sensitive enough 
to constrain the quadratic term of the dispersion relations and do not exclude 
any model in this case. Multi-messenger and multi-wavelength campaigns will 
be necessary in future to constrain this quadratic term. Combining results 
obtained with GRB and AGN detected by Fermi and ground-based detectors 
would allow considering common scheme with photon energies ranging from the 
sub-GeV range to the multi-TeV range and from sources with higher redshift 
values. Future observations with neutrino telescopes of GRB and AGN flares 
would further extend the energy range to extreme, multi-hundred TeV domain. 
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However, at present, the multi-messenger and multi-wavelength procedures are 
still to be developed. 

Concerning the increase of the experimental potential in observations of 
the AGN, a key step towards a better understanding of the emission and the 
propagation effects will be achieved with a new generation of instruments, which 
will allow the observation of a larger number of AGN flares. HESS-II, with an 
energy threshold of ^30 GeV and later the large Cherenkov array CTA0 will 
increase dramatically the sensitivity to energy-dependent time-lags. A rough 
estimation obtained by extrapolating results from the PKS 2155-304 flare to the 
extended domain in energy from a few GeV to a hundred TeV, and assuming a 
gain in sensitivity of a factor of ten, would rise the limit on the quadratic term 
by two orders of magnitude. This type of results may be further confronted with 
present and future theoretical models. It may also be underlined here that this 
will open a new unexplored domain for Lorentz Invariance Violation searches. 
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